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Abstract. A new simple Lagrangian method with favorable stabiHty and efficiency properties for computing 
general plane curve evolutions is presented. The method is based on the flowing finite volume discretization of the 
intrinsic partial differential equation for updating the position vector of evolving family of plane curves. A curve can 
be evolved in the normal direction by a combination of fourth order terms related to the intrinsic Laplacian of the 
curvature, second order terms related to the curvature, first order terms related to anisotropy and by a given external 
velocity field. The evolution is numerically stabilized by an asymptotically uniform tangential redistribution of grid 
points yielding the first order intrinsic advective terms in the governing system of equations. By using a semi-implicit 
in time discretization it can be numerically approximated by a solution to linear penta-diagonal systems of equations 
(in presence of the fourth order terms) or tri-diagonal systems (in the case of the second order terms). Various 
numerical experiments of plane curve evolutions, including, in particular, nonlinear, anisotropic and regularized 
backward curvature flows, surface diffusion and Willmore flows, are presented and discussed. 
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1. Introduction. The main purpose of this paper is to propose a simple, fast and stable 
Lagrangian method for computing evolution of closed smooth embedded plane curves driven 
by a normal velocity of the form /3(9^fc, k, x) which may depend on the intrinsic Laplacian 
d^k of the curvature k, on the curvature k itself, on the tangential angle v and the curve 
position vector x. We shall restrict our attention to the following form of the normal velocity: 

(1.1) P^-5dlk + b{k,v) + F{x). 

Here (5 > is a constant and 6 = 5(fc, v) is a smooth function satisfying 6(0, .) = 0. If (5 > 
then there is no constraint on the monotonicity of h with respect to k. On the other hand, if 
i5 = 0, then we shall assume the function 6 is strictly increasing with respect to the curvature. 

There are many interesting evolutionary models in various applied fields of science, tech- 
nology and engineering that can be described by geometric equation ( 11.11 ). For example, 
putting (5 = we obtain the normal velocity (3 = b{k, i/) + F{x) representing the well-known 
anisotropic mean curvature flow arising in the motion of material interfaces during solidifica- 
tion and in affine invariant shape analysis (see e.g. |[Tl[l5][l6l|9l[l0l[32ll22l|23l). The term 
F{x) represents an external driving force, like e.g., a given velocity field projected to the 
normal vector of a curve, or any other constant or scalar function depending on the current 
curve position. It drives the curve in the inner (if F{x) > 0) or outer (if F{x) < 0) normal 
direction. In the case S = 1 two well-known examples arise when studying a motion of the 
so-called elastic curves. It is the surface diffusion, in the case 6 = 0, and the Willmore flow 
where b = —^k^. The surface diffusion is often used in computational fluid dynamics and 
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material sciences, where the encompassing area of interface should be preserved. The case 
when b = —^k^ arises from the model of the Euler-Bernoulli elastic rod - an important prob- 
lem in structural mechanics ||6l [12] [TT| . The evolutionary models having the normal velocity 
of the form (II. lb are often adopted in image segmentation where elastic and geodesic curves 
are used in order to find image objects in an automatic way l|T9l l7l l20ll27ll29l . By our method 
we are able to handle a regularized backward mean curvature flow in which & is a decreasing 
function of the curvature k like e.g. b{k, i/) — ~k. We regularize the backward mean curva- 
ture flow by adding a small fourth order diffusion term < (5 <C 1. To our knowledge, first 
experiments of this kind are presented in this paper. 

The main idea of our approach is based on accompanying the geometric equation (II. lb by 
a stabilizing tangential velocity and in rewriting it into a form of intrinsic an partial differential 
equation (PDE) for the curve position vector The resulting PDE contains fourth, second 
and first order spatial differential terms that are approximated by means of the flowing finite 
volume method ||25l . For time discretization we follow semi-implicit approach leading to a 
solution to a linear system of equations at each time level. That can be done efficiently, and, 
due to tangential stabilization, we hope that this direct Lagrangian method can be considered 
as an efficient counterpart to the well-known level-set based methods for description of the 
curve evolution discussed in 1311 [33ll30l [TTl [131 [141 . 

The paper is organized as follows: in the next section we derive the intrinsic PDE for de- 
scription of a family of plane curves. Then we present our numerical approximation scheme. 
Finally, we discuss various numerical experiments showing applicability of our approach. 

2. Governing equations. An immersed regular plane curve F can be parameterized by 
a smooth function a; : 5^ — > M^, i.e. F = {x{u)^u G S*^} having a strictly positive local 
length term g = \dux\ > 0. Taking into account periodic boundary conditions at m = 0, 1 
we can identify the circle with the interval [0, 1]. The unit arc-length parameterization is 
denoted by s. Clearly, ds — g du. For the tangent vector we have T ^ dgX and we can choose 
the unit inward normal vector N such that det(r, N) ~ 1. The tangential angle i/ ~ arg(r), 
T = (cos ly, sinz^)^ and N = {— sin i/, cosz^)^. 

Let a regular smooth initial curve F° = Image(x'') = {x'^{u),u S [0, 1]} be given. 
We shall represent a family of plane curves F* — Image (x(.,t)) = {x{u,t),u £ [0,1]} 
that evolves according to the geometric equation dl.lb by its position vector x satisfying the 
following equation: 

(2.1) dtX = (3N + af. 

It is well-known that the presence of any tangential velocity functional a in (12.1b has no 
impact on the shape of evolving curves. However, it may help to redistribute points along the 
curve. As a consequence, it can significantly stabilize numerical computations. The reader 
is referred to papers |[I3|2l]|24l|25l|26l|27l|28l|8l|2l[3|5][35l|^ detailed discussion on 
how a suitable tangential stabilization can prevent a numerical solution from forming various 
undesired singularities. We will specify our choice of a tangential velocity a later. Since 

d^x^ d^f = a?(fciV) = d^.kN + 2dskdsN + kdlN = dlkN - 2{dsk)kf - kds{kf) 
= d^kN - 3k{dsk)f - k^d.f = d^kN - ^ds{k^)dsX - k^d^x, 

we have 



(2.2) 



{~dlk)N = -d^x - fc2a> - ^d,ik^)d,x. 



3 



Let us define the following auxiliary functions: 

(2.3) 4>{k,i^) = -Sk^ +c{k,iy), c{k,v) =b{k,v)/k 
and 

(2.4) v{k, v) = ^5 + d,<i>{k, v)-a. 

Since the function h is assumed to be smooth and 6(0, i/) = the function c(fc, v) = 
is smooth as well. Using the Frenet formulae we obtain 6(fc, v)N = c{k, v)&lx 
and v)&lx = ds{(j){k, v)dsx) — ds4>{k, v) dgX. Hence, for the normal velocity (3 of the 
form (II. lb we end up with the following higher-order intrinsic PDE for the position vector 

X = x{s, t): 

(2.5) dtx + w(fc, v)dsx = (5(-(9^.t) + ^^(^(fc, 1^)9^0;) + F{x)N{v). 

Numerical approximation of a solution x to the above PDE forms the basis of our direct 
Lagrangian approach. 

It is known (see e.g. i25 |) that a family of plane curves F* — Imagc(a;(., <)), t G [0, T), 
that evolves according to ( 12. It can be also represented by a solution to the following system 
of intrinsic parabolic-ordinary differential equations: 

(2.6) dtk = dlP + adsk + k'^l3, 

(2.7) dtv = dsf3 + ak, 

(2.8) dtg = -gkl3 + gd^a. 

In ||26ll27ll28l the system ( |2.6l l-( |278l ) was solved numerically for the curvature k, tangent angle 
V and the local length g. Knowing these quantities one can reconstruct the curve evolution. 
Asymptotically uniform redistributions for the second order flows with driving force (see 
ll26l ) and for the fourth order flows (see 1281 ) were also proposed. 

In this paper we follow a different approach when compared to 1261 |28l . It is much faster 
and simpler from computational point of view. We do not solve the system ( |2.6l )-( |278l ), but the 
position vector equation ( 12.51 ) is directly numerically discretized. On the other hand, from the 
analytical point of view, the system ( I2.6l )-( l2.8l l describes evolution of useful geometric quanti- 
ties that can be utilized in designing proper tangential velocities for stabilization of numerical 
computations. For example, in the case of convex curves, it is sufficient to solve just equation 
( 12.61 ) on the fixed parameter interval given by the range of v. Then no grid point redistribution 
is necessary l22l|23l . Equation ( 12.71 ) was used recently in designing the tangential velocities 
corresponding to the so-called crystalline curvature flow, where dtv = (c.f. l35ll34l ). In 
l26l the authors proposed and analyzed a new type of tangential redistribution referred to as 
the asymptotically uniform grid points redistribution. Here we follow a similar idea but adopt 
it in a rather different way. More precisely, in order to compute the corresponding tangential 
velocity a (depending on the curvature,tangential angle and local length) we use information 
from a solution to ( |2.5l l only, i.e. we do not solve the system ( |2.6t -( |2.8t . 

We will see that such a simple and straightforward approach yields a stable, fast and pre- 
cise solution to our problem. The corresponding numerical scheme is efficiently stabilized by 
an appropriate choice of the tangential velocity term a in ( I2.4l l. Let us note that other known 
tangential velocities, like the one preserving relative local length 118112111251 . locally diffu- 
sive redistribution |[8||27l|28|, crystalline curvature redistribution l35ll or curvature adjusted 
redistribution l34l |4l can be also incorporated straightforwardly to the numerical scheme by 
corresponding change of the term a entering equation ( I2.4l i. 
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Let us briefly outline the basic idea behind the asymptotically uniform grid points redis- 
tribution. Let us denote by Lf = Jp^ ds = g{u, t) du the length of an evolving curve F*. 
Integrating equation ( |2.81 l along the curve and taking into account periodicity of a at u = 0, 1 
we obtain 

(2.9) ^ + {kl3)TL = 0, 

where {kf3)r — x /r ^/^ denotes the curve average of the quantity kf3 over a curve T. 

When designing the tangential velocity a, it is worth to study the time evolution of the 
quantity g/L. From a numerical point of view it corresponds to the local grid point distances 
divided by the averaged local distance L /n, where n is the number of grid points. We refer 
the reader to content of the next section for details. In that sense, it represents a deviation 
of local grid point distances from being uniformly distributed. Moreover, if we define the 
quantity 9 = \n{g/L) and we take into account equations ( 12.81 ) and ( 12.91 ) we conclude 

(2.10) dtO + kf3 ^ {kp)r ^ d,a . 

From ( I2.IOI 1 we can observe that, by an appropriate choice of dga, we can control behavior 
of 9 and, subsequently, of the ratio g/L also. Our choice of a is based on the following 
particular setup (see 1261 ) 



(2.11) dsa ^ k(3 ~ {k(3)r + {e~^ ~ I) oj{t) 

where uj S Ljg^{[0, Tmax)) and Tmax is the maximal existence time of evolving curve. 
The most simple choice ui{t) = yields dt9 = 0. Hence 

9±A^9±2^ forany.E^\^e[0,T_), 

and we obtain the so-called tangential redistribution preserving the relative local length (cf. 
|[T8l|2l]|25l). On the other hand, assuming 



(2.12) / u;{T)dT = +00 

Jo 

then, by inserting (|2.11| i into ( |2.1Q| i and solving the corresponding ordinary differential equa- 
tion dt9 = (e~^ — 1) LL!{t), we obtain 9{u, — ^ as t Tmax and hence 

^1 as t ^ T„iax uniformly w.r. to it G [0, 1]. 

Lt 

It means that redistribution of grid points along a curve becomes asymptotically uniform as t 
approaches the maximal time of existence Tmax- In the case when the family {F*, t G [0, T)} 
shrinks to a point as t ^ Tmax, in order to fulfill (I2.12l i. one can choose Lo{t) = K2(fc/3}r* 
where K2 > is a positive constant. By ( I2.9l l we have w(t) dr = ~K2 /q InL'^dr = 
«;2(lnL° — In Lt) — > +cx) as t — > Tmax because limt_^T,„„^ Lt — 0. On the other hand, 
if the length Lt is always away from zero and Tmax ~ +00 one can take uj{t) ~ ki, where 
Ki > is a positive constant. Summarizing, a suitable choice of the tangential velocity 
functional a that helps to redistribute grid points uniformly along the evolved curve (in any 
case of shrinking, expanding or reaching an equilibrium curve shape) is given by a solution 
to the equation 



(2.13) dsa = kf3- {k(3)r + {L/g~l)uj, = kj + K2(fc/3)r, a{0,t) = 0, 
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where ki,K2 > 0, ki + K2 > 0, are given constants. In all computations to follow in Section 
4, we choose K2 = which is sufficient when the computations end up before an eventual 
extinction of a curve. By the boundary condition imposed on a{0, t) we prescribe the motion 
of the point x(0, t) in the normal direction only. Therefore, the tangential velocity term a is 
determined uniquely. 

It is worth to note that the speed of relaxation of the quantity 6 = ln(g/L) is controlled by 
the constant ki > 0. We typically take ki = 0(1) in all our computations. More precisely, in 
experiments presented in section 4 we take ki e [3, 10]. This is due to the fact that the speed 
of the tangential velocity a defined as in (12.13b is increasing with respect to the parameter 
Ki. Therefore choosing considerably larger values of ki would make the governing equation 
(12.51 1 strongly convective dominant which may yield a necessity of small time steps in our 
numerical scheme. 

3. Numerical approximation scheme. A family of plane curves that evolves according 
to dl.lt is numerically represented by a family of "flowing" discrete plane points where 
the index i = 1, n, denotes space discretization and the index j = 0, m, stands for 
a discrete time stepping. Assuming a uniform division of the time interval [0, T] with a 
time step t = and a uniform division of the fixed parameterization interval [0, 1] with 
a step h ~ 1/n, a discrete point x-f corresponds to x{ih,jT). Due to periodic boundary 
conditions and smoothness of the evolved curve we have used the additional points defined 
by x-'_^ — x^j^i, Xq — x^, xl^^i = x{, a;^_|_2 = The tangential velocity of a flowing 
node x-f is denoted by aj . 

The system of difference equations corresponding to (12.13b and ( 12.5b will be constructed 
at each discrete time step jr by using the flowing finite volume method proposed in |25l . First 
we solve ( 12. 13b for the tangential velocities a-f and then equation ( 12.51 ) for the position vectors 
xl,i = — 1, 71 + 2. Remaining quantities involved in (12. 13b and ( 12. 5b . as the curvature, 
tangential angle, local and total length of the curve are computed from the curve position 
vector x^^^ from the previous time step j — I. 

In order to build our numerical scheme we construct the so-called flowing finite volume 
[xl_i,xj] and also corresponding flowing dual volumes [i^, where = {xl_i+xl)/2. 
Then the approximate local lengths of flowing finite volumes , curvatures fc^ and tangential 
angles i/j are given by piecewise constant values in the flowing finite volumes. Similarly, 
the other quantities aj, xl and approximate lengths of dual volumes ql are considered as 
piecewise constant in flowing dual volumes. 

Let a discrete representation of the evolving curve be given at time level j — 1 by discrete 
points xi^"^ ,i = — 1, ?7 + 2. At every new time level j we compute 



i = 0, ...,n + 2, g; 




i+l 



) 



kl = — ^sgn(det(i?i_i,i?j+i))arccos 
2r'^ 



Rj+i-Rj-i 

i i 



) 



, j — 1 , . . . , n, /cq — /c^ , A:^ 



'II' n+l 



1 I 



vl = arccos(i?Oi/rQ), if Rq^ > , i'q = 27r — arccos(i?Oi /rg), if i?02 < , 




,i = 1, ...,n, 
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n ^ n 

1=1 1=1 



In order to compute the tangential velocity a we integrate ( 12.131 ) over the time varying 
flowing finite volume [xi-i, Xi], 



dsads= / kf3- {k/3)r + {L/g-l)ujds. 



Hereafter we use the notation J tp ds for integral of the quantity ip over the curve arc 

xt-i^i- Hence at any time level t we have the relation for approximation of the difference 

ai — a^-i 

ai - ai-i « ri{k,(3i - (fc/?)r) + [hL -r^)uj . 

Taking into account discrete time stepping in the previous relation we obtain the following 
expression for up-dated values of the tangential velocity at the j-th time level: 

(3.1) c4 =o4_^+rl{kf(3i - B^) + ihL^ -rl)uj, i^l,...,n, a^^O. 

After computing tangential velocities we update diffusion and advection terms using 

4'1 = -SiHf + c{kf,iyi), I = 1, n + 1, 



vi = — -. + aj, 1 = 1, ...,n. 

^ qi qi 

In order to compute the curve position vector x^ at a new time level we integrate ( 12.5b over 
the time varying flowing dual volume [xi, aii+i], 



dtx + v{k,iy)dsxds= / -5dfx + ds{(f>{k,iy)dsx) + F{x)N{iy) ds . 



Notice that the advective term v{k, v) contains the tangential velocity a and partial derivatives 
of functions of k and v. They can be approximated by constant difference quotients in flowing 
dual volumes. Therefore we can write 

(3.2)g,^ + - i,) = [Sdlx + 0(fc, v)dsxf:*' + qlF{xi)N {{v, + v,+{)/2) . 

dt 

Denote by d^Xi and d^xi the approximation of the k-\h arc-length derivative d^x at points Xi 
and Xi, respectively. The differences at points xi, Xi^i in the above boundary integral terms 
are then naturally approximated by 



4>{k, v)dsXi+i - (j>{k, v)dsXi « 



'^^ J 



7 



The third order terms appearing in the boundary integral in ( |3.2| | can be approximated as 
follows: 

1 / dsXi+2 - dsXi+i dsXi+i - dsXi\ 1 / dgXi+i - dgXi dgXi - d^Xi^i 



i+i \ qi+i qi / n \ g^.i 



3 3 3 J 3 3 j \ 3 3 3 3 3 3 

i 3 3 3 \ / 3 : 3 \ 



1777 Till 177 ? 7 7 7 / 

The left hand side of ( 13.2b is approximated by means of backward time differences and by a 
central finite difference approximation of the advective term (up-wind technique can be also 
easily incorporated). We obtain 

at T I 

Taking into account information from the previous time step j — 1 in the driving term FN ^ 
multiplying the third order terms approximation by — 5 and putting all unknowns xj, I = 
i ~ 2, . . . ,i + 2 with their coefficients to the left hand side, we obtain a linear system of 
equations for updating the discrete position vector x^ of the evolved curve T^'^: 



for i — 1, n, subject to periodic boundary conditions x-'_i — xl^_i, Xq — xl^, x-' 



n+l 



, xl^^2 = ^2 where 



3 3 3 3 3 3 



V: 




^ ^xf' + qi Fixj-^)N 

The previous system is penta-diagonal if (5 > and tri-diagonal in the case S = 0. In the 
latter case when (5 = 0, the monotonicity assumption on the function b guarantees the strict 
diagonal dominance of the tri-diagonal system matrix. In both cases, it is solved efficiently by 
means of the Gauss-Seidel iterative method or by its well-known successive over relaxation 
version SOR. We start the iterates from the previous time step vector x^~^ and, in practice, 
when using our asymptotically uniform tangential redistribution (AUTR), there are just few 
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b) 



Fig. 4.1. Affine invariant shrinking evolution of an initial ellipse (b = k^/^,& = 0, F = 0); a) with 
asymptotically uniform redistribution (k\ = 3j,' b) without redistribution. Numerical parameters n = 100, t = 
0.001. Time steps j = 0, 200, ■ • • , 1400 are plotted using lines and grid points, while time steps j = 100, 300, ■ • ■ , 
up to j = 1500 a); and j = 1100 b) are plotted using lines. 

number of SOR iterations needed in order to achieve the solver accuracy goal. The iterative 
process is stopped when a difference of subsequent iterates in the maximum norm is less 
than the prescribed tolerance, e.g. TOL = lO^^*'. Based on our practical experience, it 
should be also noted that the number of SOR iterations (corresponding to well-conditioning 
of an iteration matrix) is strongly lowered by such a proper choice of the tangential velocity. 
Consequently, it speeds up computations significantly making thus our numerical scheme 
fast and efficient. Computational times are reported in next section. In some examples of 
evolution of curves with high variation in the curvature or in checking the experimental order 
of convergence (see e.g. Fig. 14.71 or Tables WA\ and |4~2] ) we typically use small time steps 
T K, h? or even t « /i"' (in some nontrivial cases of the fourth order flows). Nevertheless, 
taking larger time steps t leads to satisfactory numerical results as it can be seen from other 
experiments presented in this paper We observed neither occurrences of accumulation of 
grid points nor spurious numerical oscillations. 

4. Discussion on numerical experiments. In the first numerical experiments we show 
a stabilizing effect of our scheme due to asymptotically uniform tangential redistribution 
(AUTR) for the case of selfsimilar affine invariant shrinking evolution of an initial ellipse with 
half-axes ratio 3:1 (Fig. 14. lb ). When the grid points are moving only in the normal direction, 
the numerical computation collapses soon because of merging of grid points and spurious 
swallow-tails creation in the left and right end of the ellipse (Fig. 14.1b ). In Fig. 14.21 we 
present an anisotropic curve shortening evolution of an ellipse computed again with help of 
the asymptotically uniform tangential redistribution. In all numerical experiments parameters 
of computations are shown in figure captions. 

Again starting from an initial ellipse with half-axes ratio 1:3 we numerically compute 
its evolution by the surface diffusion. Fig. 14.3b shows the result with and Fig. 14.3b without 
asymptotically uniform tangential redistribution (AUTR). In Tables l4n and l4.2l we show how 
precisely the encompased area is preserved (which is one of the analytical properties of sur- 
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Fig. 4.2. Anisotropic shrinking evolution of an initial ellipse with AJJTR; 6 = (1 — 0.9 cos(4i/ — tt)) k,5 = 
0, F = 0. Numerical andAUTR parameters: n = 100, r = 0.001, ki = 3. Time steps j = 0, 200, ■ • ■ , 1400 are 
plotted using lines and grid points, while time steps j = 100, 300, • ■ ■ , 1500 are plotted using lines. 




Fig. 4.3. Surface diffusion flow (<5=l,fe = 0,-F = 0j of an initial ellipse with a) and without tangential 
redistribution b). Numerical and redistribution parameters: n = 100, t = 0.001, ki = 10. Time steps j = 
0, 400, ■ ■ ■ , 2000 are plotted using lines and grid points, while time steps j = 200, 600, ■ • • , 1800 are plotted 
using lines. 

Table 4.1 

An ellipse evolving by the surface diffusion using AUTR, same parameters as in Fig. \4.3U - We report errors in 
area evolution, experimental order of convergence (EOC) in this quantity and computational time (CPU) for refined 
discretization parameters. 



n 


T 


# of steps 


area error 


EOC 


CPU (sec) 


25 


0.016 


125 


0.0774 




0.02 


50 


0.004 


500 


0.0202 


1.93 


0.43 


100 


0.001 


2000 


0.0052 


1.89 


2.94 


200 


0.00025 


8000 


0.0014 


1.91 


40.21 



Table 4.2 

An ellipse evolving by the surface diffusion without tangential redistribution, same parameters as in Fig. \4.3b . 
Same quantities as above are reported. 



n 


T 


# of steps 


area error 


EOC 


CPU (sec) 


25 


0.016 


125 


0.2922 




0.11 


50 


0.004 


500 


0.1056 


1.46 


4.66 


100 


0.001 


2000 


0.0320 


1.72 


198 


200 


0.00025 


8000 


0.0155 


1.04 


1741 
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face diffusion flow with (5 = 1 and 6 = 0) during the evolution. The main observation here 
is the fact that when using AUTR we obtain significantly lower error and the method with 
AUTR has approximately the second order of accuracy which is not the case for computa- 
tions without AUTR. In these experiments we used coupling r « /i^, cf. ll9l [T0l[mi27ll28l . 
the exact area at any time moment is = Stt and the area error is computed as 

(4.1) ||CI!= (E(^'-^«)'^ 

where = | dct(a;:^ , — gives the encompassed area of the polygonal com- 

puted curve at the th time step. 

In Tables l4n and l4.2l we also report computational times achieved on a standard 2.2GHz 
laptop. They justify why the method with AUTR is refereed to as a fast method. For standard 
curve resolutions, with 100 or 200 grid points, we get that one time step takes 0.0014 respec- 
tively 0.0050 second. Such fast CPU times are obtained due to tangential redistribution which 
not only stabilized the computations but also improve cyclic penta-diagonal iteration matrix 
properties in the SOR iterative method (the relaxation parameter was set to 1.6). Such CPU 
times are obtained also in further experiments presented in this section. Of course, one has 
to multiply them by the number of time steps which may lead to overall long computations 
when computing long time behaviors of curve evolutions shown e.g. in Fig. 14. 81 

The next set of experiments is focused on the backward mean curvature flow with ex- 
panding constant force. The flow is regularized by various strengths of the Willmore flow. 
Starting with an initial ellipse with half-axes ratio 3:1 we can observe that in the case of 
a strong regularization the backward mean curvature flow is dominated by the elastic relax- 
ation due to the Willmore surface energy (see Fig. 14.4b ). On the other hand, if we decrease the 
fourth order regularization by taking smaller values of the parameter S > 0, the nonconvex 
parts formed during evolution are attenuated (see Fig. 14. 4b ) and even curve selfintersections 
may occur (see Figs. 14.51 and |43] l as it can be expected in the backward in time diffusion 
process. 

In Fig. |4.5| we also show an important role of AUTR in such nontrivial experiments. The 
asymptotically uniform tangential redistribution keeps very good curve resolution even in 
cases of several subsequent curve selfintersections. On the other hand, if we consider for ex- 
ample the well-known redistribution preserving the relative local length ifTSl 12111251 . obtained 
by taking = in ( |2.11| i and (13. It then this method is not capable to handle this situations 
properly. If we consider 6 = 0.01 the backward diffusion effects are strongly dominating. In 
Fig. 14.61 we can observe very fast shape "enrichment" soon after even a nonconvexification 
of the evolving curve. It is worth to note that such experiments would be impossible without 
incorporating the tangential stabilization into the direct Lagrangian computational approach. 

The last set of experiments is devoted to evolution of an initial spiral given by 

Xi{u)=acosb, X2{u)~asuib, 

a = 0.5 e"^"5 ^in{27ru) _ Q cos(27rit), b = 10 arctan(l + 0.5 sin(27rw)) , 

plotted in Fig. 14. SB . Again the presence of the tangential redistribution (AUTR) has a stabiliz- 
ing effect on all numerical computations. Without redistribution a parametric approach col- 
lapses soon. The evolution by the mean curvature and surface diffusion is shown in Fig. 14.7b 
and Fig. 14.7b . The backward curve diffusion regularized by a different strength of the Will- 
more flow is presented in Fig. l4.8b -c. 
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a) b) 



Fig. 4.4. Backward mean cun'ature flow with negative external force regularized by the Willmore surface 
energy (b = —k — ^5k^,F = —1). a) a strong regularization effect with 5 = 1; h) a weak regularization 
& = 0.1. Numerical and AUTR parameters: n = 100, r = 0.001, ki = 10. In subfigure a) the time steps 
j = 0, 400, ■ ■ ■ , 2000 are plotted using lines and grid points, while time steps j = 200, 600, • ■ ■ , 1800 are plotted 
only with lines. In subfigure b) the time steps j = 0, 200, • ■ • , 1800 are plotted. 




Fig. 4.5. Comparison of continuation of the evolution from Fieure \4M b} computed with AUTR (upper row) 
and with redistribution preserving relative local length (bottom row), respectively. In both cases we show time steps 
j = 2000 (left) and j = 2200 (right). 




Fig. 4.6. Regularized backward curvature driven flow with b = ~k — ^Sk^,F = —1 and even weaker 
regularization with 5 = 0.01. Redistribution parameter: ki = 10 and numerical parameters: n = 400, t = 
0.0001, are used. Time steps j = 0, 1000, 2000, 3000 are plotted. 
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a) b) 



Fig. 4.7. a) mean curvature flow with b(k, u) = k, S = 0, F = 0; b) surface diffusion flow with b = 0, 5 = 
l.F = of an initial spiral. Numerical and AUTR parameters: n = 100, ki = 10 and t = 10~^ a) and 
T = 10-10 h). Time steps j G {0,1,2,3,4,5,6,6.5} X a) and j G {0,1,5,10,20,40,60,80,100} X W 
bj are plotted. 

5. Conclusions. A new direct Lagrangian method stabilized by a suitable tangential re- 
distribution has been presented for the case of general plane curve evolution models. An 
evolved curve is driven in the normal direction by a combination of the fourth order terms re- 
lated to the intrinsic Laplacian of curvature, second order terms related to the curvature, first 
order terms related to anisotropy and tangential redistribution and by a given external veloc- 
ity field. We showed how a proper choice of a tangential velocity can stabilize and speed up 
computations. Nontrivial numerical experiments justified applicability and numerical stabil- 
ity of the approximation scheme in the anisotropic mean curvature flow, surface diffusion and 
the Willmore flow and even in the case of backward curve diffusion slightly regularized by 
the fourth order terms. Further study of the scheme from the numerical analysis point of view 
as well as derivation of analytical quantities to which the numerical results can be compared 
in the above mentioned models will be the objective of our future research. 
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Fig. 4.8. Bachvard mean cutMature flow regularized by the Willmore flow (b(k, u) = —k — ^&k^ , F = 0) 
with diflerent strengths 5 = I a), 5 = 0.1 h), 5 = 0.01 c) of the same initial spiral d). Numerical and AUTR 
parameters: n = 200, t = 10"^^, ki = 100. Time steps: j = 0, lO'^ , 10*^, 10^, 10*, 4 10*, lO'' are p/ofterf. 
Initial and the last step plotted with grid points. 
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